function tfc_powered_descent()
    %========================
    % 1️⃣ 设置初始条件
    %========================
    Tmax = 13258.18;   % 最大推力 [N]
    Tmin = 4971.81;    % 最小推力 [N]
    Isp = 225;         % 比冲 [s]
    g0 = 9.807;        % 地球重力加速度 [m/s^2]
    ag = [0; 0; -3.7114];  % 火星重力 [m/s^2]
    m0 = 2000;         % 初始质量 [kg]

    % 初始与目标状态
    r0 = [2000; 0; 1000];  % 初始位置 [m]
    v0 = [-50; 0; -20];    % 初始速度 [m/s]
    rf = [0; 0; 0];        % 目标位置 [m]
    vf = [0; 0; 0];        % 目标速度 [m/s]

    % 设定时间参数
    % 设定离散点数
    discretize_points = 100;
    tf_guess = 50; % 着陆时间
    t = linspace(0, tf_guess, discretize_points);

    %========================
    % 2️⃣ 计算 TFC 约束表达式
    %========================
    N = 10;  % Chebyshev 多项式阶数
    [TFC_mat, basis] = tfc_basis_functions(N, t, tf_guess);

    %========================
    % 3️⃣ 最小二乘法优化轨迹
    %========================
    coeff = lsqnonlin(@(c) tfc_residuals(c, TFC_mat, basis, t, r0, v0, rf, vf, m0, Tmax, Tmin, Isp, ag), ...
                      zeros(3 * N, 1));

    %========================
    % 4️⃣ 计算最优轨迹并绘制
    %========================
    r_opt = TFC_mat * coeff(1:N);
    v_opt = TFC_mat * coeff(N+1:2*N);
    
    figure;
    plot3(r_opt(:,1), r_opt(:,2), r_opt(:,3), 'r', 'LineWidth', 2); hold on;
    plot3(rf(1), rf(2), rf(3), 'bo', 'MarkerSize', 8, 'MarkerFaceColor', 'b');
    xlabel('X [m]'); ylabel('Y [m]'); zlabel('Z [m]');
    title('TFC Powered Descent Trajectory');
    grid on;
end

%========================
% 🚀 2.1 计算 TFC 约束表达式
%========================
function [TFC_mat, basis] = tfc_basis_functions(N, t, tf)
    basis = chebyshev_basis(N, t / tf);
    TFC_mat = [basis, basis .* t, basis .* t.^2];
end

%========================
% 🚀 2.2 计算 Chebyshev 多项式基函数
%========================
function basis = chebyshev_basis(N, t)
    basis = zeros(length(t), N);
    for k = 1:N
        basis(:,k) = cos((k-1) * acos(2*t - 1));
    end
end

%========================
% 🚀 2.3 计算 TFC 残差函数
%========================
function res = tfc_residuals(coeff, TFC_mat, basis, t, r0, v0, rf, vf, m0, Tmax, Tmin, Isp, ag)
    r_coeff = coeff(1:length(basis));
    v_coeff = coeff(length(basis)+1:2*length(basis));

    % 计算位置和速度
    r = TFC_mat * r_coeff;
    v = TFC_mat * v_coeff;

    % 计算动力学残差
    r_dot = v;
    v_dot = ag + optimal_thrust_tfc(t, Tmax, Tmin) / m0 .* (-v ./ vecnorm(v, 2, 2));

    res = [r_dot - v; v_dot - ag];
end

%========================
% 🚀 2.4 最优推力策略 (Max-Min-Max)
%========================
function T = optimal_thrust_tfc(t, Tmax, Tmin)
    T = Tmax * (t < 15) + Tmin * (t >= 15 & t < 35) + Tmax * (t >= 35);
end
